
import openseespy.opensees as ops

ops.wipe()

ops.model('basic', '-ndm', 3, '-ndf', 6)

###################################
# Define Material
###################################

# Define PSUMAT and convert it to plane stress material
ops.nDMaterial('PlaneStressUserMaterial', 1, 40, 7, 20.7e6, 2.07e6, -4.14e6, -0.002, -0.01, 0.001, 0.3)
ops.nDMaterial('PlateFromPlaneStress', 4, 1, 1.25e10)

# Define material for rebar
ops.uniaxialMaterial('Steel02', 7, 379e6, 202.7e9, 0.01, 18.5, 0.925, 0.15)
ops.uniaxialMaterial('Steel02', 8, 392e6, 200.6e9, 0.01, 18.5, 0.925, 0.15)

# Convert rebar material to plane stress/plate rebar
# Angle 0 is for vertical rebar and 90 is for horizontal rebar
ops.nDMaterial('PlateRebar', 9, 7, 90.0)
ops.nDMaterial('PlateRebar', 10, 8, 90.0)
ops.nDMaterial('PlateRebar', 11, 8, 0.0)

# Define LayeredShell sections. Section 1 is used for the special boundary elements and section 2 is used for
# the unconfined interior wall portion
ops.section('LayeredShell', 1, 10, 4, 0.0125, 11, 0.0002403, 11, 0.0003676, 4, 0.024696, 4, 0.024696, 4,
            0.024696, 4,
            0.024696, 11, 0.0003676, 11, 0.0002403, 4, 0.0125)
ops.section('LayeredShell', 2, 8, 4, 0.0125, 11, 0.0002403, 10, 0.0002356, 4, 0.0495241, 4, 0.0495241, 10,
            0.0002356,
            11, 0.0002403, 4, 0.0125)

# ##################
# NODES
# ##################
# define nodes
for i in range(11):
    ops.node(16 * i + 1, 0, 0, 0.2 * i)
    ops.node(16 * i + 2, 0.2, 0, 0.2 * i)
    ops.node(16 * i + 3, 0.5, 0, 0.2 * i)
    ops.node(16 * i + 4, 0.8, 0, 0.2 * i)
    ops.node(16 * i + 5, 1.0, 0, 0.2 * i)

    ops.node(16 * i + 6, 0, 1.0, 0.2 * i)
    ops.node(16 * i + 7, 0.2, 1, 0.2 * i)
    ops.node(16 * i + 8, 0.5, 1, 0.2 * i)
    ops.node(16 * i + 9, 0.8, 1, 0.2 * i)
    ops.node(16 * i + 10, 1.0, 1, 0.2 * i)

    ops.node(16 * i + 11, 0, 0.2, 0.2 * i)
    ops.node(16 * i + 12, 0, 0.5, 0.2 * i)
    ops.node(16 * i + 13, 0, 0.8, 0.2 * i)
    ops.node(16 * i + 14, 1, 0.2, 0.2 * i)
    ops.node(16 * i + 15, 1, 0.5, 0.2 * i)
    ops.node(16 * i + 16, 1, 0.8, 0.2 * i)

##########################
# ELEMENTS
##########################

ShellType = "ShellNLDKGQ"
# ShellType = "ShellMITC4"
for i in range(10):
    ops.element(ShellType, i + 1, 16 * i + 1, 16 * i + 2, 16 * (i + 1) + 2, 16 * (i + 1) + 1, 1)
    ops.element(ShellType, i + 11, 16 * i + 2, 16 * i + 3, 16 * (i + 1) + 3, 16 * (i + 1) + 2, 2)
    ops.element(ShellType, i + 21, 16 * i + 3, 16 * i + 4, 16 * (i + 1) + 4, 16 * (i + 1) + 3, 2)
    ops.element(ShellType, i + 31, 16 * i + 4, 16 * i + 5, 16 * (i + 1) + 5, 16 * (i + 1) + 4, 1)

    ops.element(ShellType, i + 41, 16 * i + 6, 16 * i + 7, 16 * (i + 1) + 7, 16 * (i + 1) + 6, 1)
    ops.element(ShellType, i + 51, 16 * i + 7, 16 * i + 8, 16 * (i + 1) + 8, 16 * (i + 1) + 7, 2)
    ops.element(ShellType, i + 61, 16 * i + 8, 16 * i + 9, 16 * (i + 1) + 9, 16 * (i + 1) + 8, 2)
    ops.element(ShellType, i + 71, 16 * i + 9, 16 * i + 10, 16 * (i + 1) + 10, 16 * (i + 1) + 9, 1)

    ops.element(ShellType, i + 81, 16 * i + 1, 16 * i + 17, 16 * i + 27, 16 * i + 11, 1)
    ops.element(ShellType, i + 91, 16 * i + 11, 16 * i + 27, 16 * i + 28, 16 * i + 12, 2)
    ops.element(ShellType, i + 101, 16 * i + 12, 16 * i + 28, 16 * i + 29, 16 * i + 13, 2)
    ops.element(ShellType, i + 111, 16 * i + 13, 16 * i + 29, 16 * i + 22, 16 * i + 6, 1)

    ops.element(ShellType, i + 121, 16 * i + 5, 16 * i + 14, 16 * i + 30, 16 * i + 21, 1)
    ops.element(ShellType, i + 131, 16 * i + 14, 16 * i + 15, 16 * i + 31, 16 * i + 30, 2)
    ops.element(ShellType, i + 141, 16 * i + 15, 16 * i + 16, 16 * i + 32, 16 * i + 31, 2)
    ops.element(ShellType, i + 151, 16 * i + 16, 16 * i + 10, 16 * i + 26, 16 * i + 32, 1)

# P-delta columns

for i in range(10):
    ops.element('truss', 500 + i + 1, 16 * i + 1, 16 * (i + 1) + 1, 223.53e-6, 7)
    ops.element('truss', 500 + i + 11, 16 * i + 2, 16 * (i + 1) + 2, 223.53e-6, 7)
    ops.element('truss', 500 + i + 21, 16 * i + 4, 16 * (i + 1) + 4, 223.53e-6, 7)
    ops.element('truss', 500 + i + 31, 16 * i + 5, 16 * (i + 1) + 5, 223.53e-6, 7)

    ops.element('truss', 500 + i + 41, 16 * i + 6, 16 * (i + 1) + 6, 223.53e-6, 7)
    ops.element('truss', 500 + i + 51, 16 * i + 7, 16 * (i + 1) + 7, 223.53e-6, 7)
    ops.element('truss', 500 + i + 61, 16 * i + 9, 16 * (i + 1) + 9, 223.53e-6, 7)
    ops.element('truss', 500 + i + 71, 16 * i + 10, 16 * (i + 1) + 10, 223.53e-6, 7)

    ops.element('truss', 500 + i + 81, 16 * i + 11, 16 * (i + 1) + 11, 223.53e-6, 7)
    ops.element('truss', 500 + i + 91, 16 * i + 13, 16 * (i + 1) + 13, 223.53e-6, 7)

    ops.element('truss', 500 + i + 101, 16 * i + 14, 16 * (i + 1) + 14, 223.53e-6, 7)
    ops.element('truss', 500 + i + 111, 16 * i + 16, 16 * (i + 1) + 16, 223.53e-6, 7)

for i in range(176):
    m = 246000 * 4 / 9.8 / 176
    ops.mass(i + 1, m, m, m, 0, 0, 0)
# Fix all bottom nodes
ops.fixZ(0, 1, 1, 1, 1, 1, 1)